rm(list = ls())
#working_folder<- "/Users/bgpopescu/Library/CloudStorage/Dropbox/covid_institutions/"
working_folder<-"//Mac/Dropbox-1/covid_paper_replication/"
setwd(working_folder)

#This is to obtain:
#-figure_2

library("readxl")
library("reshape")
library("stringi")
library("stringr")
library("dplyr")
library("tidyr")
library("stargazer")
library("corrplot")
library("Hmisc")
library("glue")
library("ggplot2")
library("glue")
library('gridExtra')
library('grid')
library('lfe')
library('broom')
library("broom")
library("ggrepel")
library("lemon")
library("haven")
library("plm")
library("lmtest")
library("multiwayvcov")
library("haven")
library("data.table")
library("sf")


##################
#Step 2 Bandwidth#
##################

lits_shape_2016 <- st_read(dsn="./data/shape_files.gdb",
                             layer="lits_shape_2016")

data_2016<-lits_shape_2016
data_2016$treat <- as.numeric( with (data_2016,ifelse((data_2016$treat==1), 1, 0)))
data_2016$dist_east_west_brd<-data_2016$east_west_distance/1000
data_2016$dist1 <- as.numeric( with (data_2016,ifelse(data_2016$treat==1, 1, -1)))
data_2016$dist2 <-data_2016$dist1*data_2016$dist_east_west_brd

data_2016$bfe1<-0
data_2016$bfe1[data_2016$east_west_NEAR_FID==4]<-1
data_2016$bfe2<-0
data_2016$bfe2[data_2016$east_west_NEAR_FID==5]<-1
data_2016$bfe3<-0
data_2016$bfe3[data_2016$east_west_NEAR_FID==6]<-1


lits_shape_2010 <- st_read(dsn="./data/shape_files.gdb",
                           layer="lits_shape_2010")

data_2010<-lits_shape_2010
data_2010$treat <- as.numeric( with (data_2010,ifelse((data_2010$treat==1), 1, 0)))
data_2010$dist_east_west_brd<-data_2010$east_west_distance/1000
data_2010$dist1 <- as.numeric( with (data_2010,ifelse(data_2010$treat==1, 1, -1)))
data_2010$dist2 <-data_2010$dist1*data_2010$dist_east_west_brd

data_2010$bfe1<-0
data_2010$bfe1[data_2010$east_west_NEAR_FID==4]<-1
data_2010$bfe2<-0
data_2010$bfe2[data_2010$east_west_NEAR_FID==5]<-1
data_2010$bfe3<-0
data_2010$bfe3[data_2010$east_west_NEAR_FID==6]<-1



###########
#LITS 2016#
###########

list2016<-read_dta("./data/2016-LiTS-III.dta", encoding = 'latin1')
Ger_2016<- list2016[ which(list2016$country=="Germany"), ]
rm(list2016)
names(Ger_2016)
Ger_2016$region_name


###########
#LITS 2010#
###########

list2010<-read_dta("./data/2010 - lits2.dta", encoding = 'latin1')
Ger_2010<- list2010[ which(list2010$country_=="Germany"), ]
rm(list2010)
names(Ger_2010)
Ger_2010$Region1


############################################################
#SPORT AND RECREATIONAL ORGANIZATIONS AND ASSOCIATIONS 2016#
############################################################

Ger_2016$q919b
unique(Ger_2016$q919b)

#GR1
Ger_2016$member_sport<-Ger_2016$q919b
Ger_2016$member_sport[Ger_2016$q919b==1] <- 1
Ger_2016$member_sport[Ger_2016$q919b==2] <- 1
Ger_2016$member_sport[Ger_2016$q919b==3] <- 0

#GR2
Ger_2016$member_sport_active_passive<-NA
Ger_2016$member_sport_active_passive[Ger_2016$q919b==1] <- 1
Ger_2016$member_sport_active_passive[Ger_2016$q919b==2] <- 0

#GR3
Ger_2016$member_sport_active_nomb<-NA
Ger_2016$member_sport_active_nomb[Ger_2016$q919b==1] <- 1
Ger_2016$member_sport_active_nomb[Ger_2016$q919b==3] <- 0

#GR4
Ger_2016$member_sport_active_all<-NA
Ger_2016$member_sport_active_all[Ger_2016$q919b==1] <- 1
Ger_2016$member_sport_active_all[Ger_2016$q919b==2] <- 0
Ger_2016$member_sport_active_all[Ger_2016$q919b==3] <- 0

#GR5
Ger_2016$member_sport_passive_active<-NA
Ger_2016$member_sport_passive_active[Ger_2016$q919b==2] <- 1
Ger_2016$member_sport_passive_active[Ger_2016$q919b==1] <- 0

#GR6
Ger_2016$member_sport_passive_nonmb<-NA
Ger_2016$member_sport_passive_nonmb[Ger_2016$q919b==2] <- 1
Ger_2016$member_sport_passive_nonmb[Ger_2016$q919b==3] <- 0

#GR7
Ger_2016$member_sport_passive_all<-NA
Ger_2016$member_sport_passive_all[Ger_2016$q919b==2] <- 1
Ger_2016$member_sport_passive_all[Ger_2016$q919b==1] <- 0
Ger_2016$member_sport_passive_all[Ger_2016$q919b==3] <- 0


summary(Ger_2016$member_sport_active)
summary(Ger_2016$member_sport_active_nomb)
summary(Ger_2016$member_sport_active_all)
summary(Ger_2016$member_sport_passive)
summary(Ger_2016$member_sport_passive_active)
summary(Ger_2016$member_sport_passive_nonmb)
summary(Ger_2016$member_sport_passive_all)


############################################################
#SPORT AND RECREATIONAL ORGANIZATIONS AND ASSOCIATIONS 2010#
############################################################

#GR1
Ger_2010$member_sport<-Ger_2010$q713b
Ger_2010$member_sport[Ger_2010$q713b==1] <- 1
Ger_2010$member_sport[Ger_2010$q713b==2] <- 1
Ger_2010$member_sport[Ger_2010$q713b==3] <- 0

#GR2
Ger_2010$member_sport_active_passive<-NA
Ger_2010$member_sport_active_passive[Ger_2010$q713b==1] <- 1
Ger_2010$member_sport_active_passive[Ger_2010$q713b==2] <- 0

#GR3
Ger_2010$member_sport_active_nomb<-NA
Ger_2010$member_sport_active_nomb[Ger_2010$q713b==1] <- 1
Ger_2010$member_sport_active_nomb[Ger_2010$q713b==3] <- 0

#GR4
Ger_2010$member_sport_active_all<-NA
Ger_2010$member_sport_active_all[Ger_2010$q713b==1] <- 1
Ger_2010$member_sport_active_all[Ger_2010$q713b==2] <- 0
Ger_2010$member_sport_active_all[Ger_2010$q713b==3] <- 0

#GR5
Ger_2010$member_sport_passive_active<-NA
Ger_2010$member_sport_passive_active[Ger_2010$q713b==2] <- 1
Ger_2010$member_sport_passive_active[Ger_2010$q713b==1] <- 0

#GR6
Ger_2010$member_sport_passive_nonmb<-NA
Ger_2010$member_sport_passive_nonmb[Ger_2010$q713b==2] <- 1
Ger_2010$member_sport_passive_nonmb[Ger_2010$q713b==3] <- 0

#GR7
Ger_2010$member_sport_passive_all<-NA
Ger_2010$member_sport_passive_all[Ger_2010$q713b==2] <- 1
Ger_2010$member_sport_passive_all[Ger_2010$q713b==1] <- 0
Ger_2010$member_sport_passive_all[Ger_2010$q713b==3] <- 0


summary(Ger_2010$member_sport)
summary(Ger_2010$member_sport_active_passive)
summary(Ger_2010$member_sport_active_nomb)
summary(Ger_2010$member_sport_active_all)
summary(Ger_2010$member_sport_passive_active)
summary(Ger_2010$member_sport_passive_nonmb)
summary(Ger_2010$member_sport_passive_all)


####################
#Law obeyance 2016#
####################
Ger_2016$law_obeyance<-Ger_2016$q417d
Ger_2016$law_obeyance[Ger_2016$law_obeyance==-97] <- NA

####################
#Law obeyance 2010#
####################
Ger_2010$law_obeyance<-Ger_2010$q316d
Ger_2010$law_obeyance[Ger_2010$law_obeyance==-97] <- NA


#####
#Age#
#####


#Age 2016
Ger_2016$age_group<-ifelse(Ger_2016$age_pr>=18 & Ger_2016$age_pr<=24,  1,
                           ifelse(Ger_2016$age_pr>=25 & Ger_2016$age_pr<=34, 2,
                                  ifelse(Ger_2016$age_pr>=35 & Ger_2016$age_pr<=44, 3,
                                         ifelse(Ger_2016$age_pr>=45 & Ger_2016$age_pr<=54, 4,
                                                ifelse(Ger_2016$age_pr>=65, 5, NA)))))
#Age 2010
Ger_2010$age_group<-ifelse(Ger_2010$q104a_1>=18 & Ger_2010$q104a_1<=24,  1,
                           ifelse(Ger_2010$q104a_1>=25 & Ger_2010$q104a_1<=34, 2,
                                  ifelse(Ger_2010$q104a_1>=35 & Ger_2010$q104a_1<=44, 3,
                                         ifelse(Ger_2010$q104a_1>=45 & Ger_2010$q104a_1<=54, 4,
                                                ifelse(Ger_2010$q104a_1>=65, 5, NA)))))

#Education 2016
Ger_2016$edu<-ifelse(Ger_2016$q109_1==1, 1,
                     ifelse(Ger_2016$q109_1==2, 2,
                            ifelse(Ger_2016$q109_1==3 | Ger_2016$q109_1==4, 3,
                                   ifelse(Ger_2016$q109_1==5, 4,
                                          ifelse(Ger_2016$q109_1==6, 5,
                                                 ifelse(Ger_2016$q109_1==7, 6, NA))))))
#Education 2010
Ger_2010$q515
Ger_2010$edu<-ifelse(Ger_2010$q515==1, 1,
                     ifelse(Ger_2010$q515==2, 2,
                            ifelse(Ger_2010$q515==3 | Ger_2010$q515==4, 3,
                                   ifelse(Ger_2010$q515==5, 4,
                                          ifelse(Ger_2010$q515==6, 5,
                                                 ifelse(Ger_2010$q515==7, 6, NA))))))


#Income
Ger_2016$income<-as.character(as.factor(Ger_2016$PRq315))
Ger_2016$income<-ifelse(Ger_2016$income=="-97" | Ger_2016$income=="-99", NA, Ger_2016$income)
unique(Ger_2016$income)
Ger_2016$income<-as.numeric(Ger_2016$income)
Ger_2016$PSU_name
Ger_2016$weight_sample
Ger_2016$weight_population

Ger_2010$income<-Ger_2010$q227
Ger_2010$income<-ifelse(Ger_2010$income=="-97" | Ger_2010$income=="-99", NA, Ger_2010$income)
Ger_2010$psu
Ger_2010$weight
Ger_2010$Region1
Ger_2010$SerialID

Ger_2010_b<-subset(Ger_2010, select = c(SerialID,
                                        Region1,
                                        member_sport,
                                        member_sport_active_passive,
                                        member_sport_active_nomb,
                                        member_sport_active_all,
                                        member_sport_passive_active,
                                        member_sport_passive_nonmb,
                                        member_sport_passive_all,
                                        law_obeyance,
                                        age_group,
                                        edu,
                                        income,
                                        psu))


Ger_2016_b<-subset(Ger_2016, select = c(ID,
                                        region_name,
                                        member_sport,
                                        member_sport_active_passive,
                                        member_sport_active_nomb,
                                        member_sport_active_all,
                                        member_sport_passive_active,
                                        member_sport_passive_nonmb,
                                        member_sport_passive_all,
                                        law_obeyance,
                                        age_group,
                                        edu,
                                        income,
                                        PSU_name))


Ger_2016_b$ID<-as.character(as.factor(Ger_2016_b$ID))
unique(duplicated(Ger_2016_b$ID))
unique(duplicated(data_2016$ID))


Ger_2016_b$ID
data_2016$ID<-as.character(data_2016$ID)

data_2016<-left_join(data_2016, Ger_2016_b, by = c("ID"="ID"))
data_2016$POINT_Y<-data_2016$latitude
data_2016$POINT_X<-data_2016$longitude
data_2016$dist_east_west_brd
data_2016$PSU_name
data_2016$region_name

data_2010$SerialID
Ger_2010_b$SerialID

data_2010<-left_join(data_2010, Ger_2010_b, by = c("SerialID"="SerialID"))
data_2010$POINT_Y<-data_2010$latitude
data_2010$POINT_X<-data_2010$longitude
data_2010$dist_east_west_brd
data_2010$PSU_name<-as.character(data_2010$psu)
data_2010$region_name<-data_2010$Region1



data_2016$wave<-2016
data_2010$wave<-2010
combined<-dplyr::bind_rows(data_2010, data_2016)




#####################
#SPORT: Law obeyance#
#####################
#GR1
law_obeyance_sport_2016<-lm(law_obeyance~member_sport+
                              age_group +income +edu+
                              factor(region_name),
                            data=data_2016, 
                            na.action=na.exclude)
summary(law_obeyance_sport_2016)
law_obeyance_sport_2016_se <- vcovHC(law_obeyance_sport_2016, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_2016_se <- sqrt(diag(law_obeyance_sport_2016_se))

#GR2
law_obeyance_sport_active_passive_2016<-lm(law_obeyance~member_sport_active_passive+
                                             age_group +income +edu+
                                             factor(region_name),
                                           data=data_2016, na.action=na.exclude)
summary(law_obeyance_sport_active_passive_2016)
law_obeyance_sport_active_passive_2016_se <- vcovHC(law_obeyance_sport_active_passive_2016, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_active_passive_2016_se <- sqrt(diag(law_obeyance_sport_active_passive_2016_se))

#GR3
law_obeyance_sport_active_nomb_2016<-lm(law_obeyance~member_sport_active_nomb+
                                          age_group +income +edu+
                                          latitude + longitude+treat+
                                          factor(region_name),
                                        data=data_2016, 
                                        na.action=na.exclude)
summary(law_obeyance_sport_active_nomb_2016)
law_obeyance_sport_active_nomb_2016_se <- vcovHC(law_obeyance_sport_active_nomb_2016, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_active_nomb_2016_se <- sqrt(diag(law_obeyance_sport_active_nomb_2016_se))

#GR4
law_obeyance_sport_active_all_2016<-lm(law_obeyance~member_sport_active_all+
                                         age_group +income +edu+
                                         factor(region_name),
                                       data=data_2016, na.action=na.exclude)
summary(law_obeyance_sport_active_all_2016)
law_obeyance_sport_active_all_2016_se <- vcovHC(law_obeyance_sport_active_all_2016, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_active_all_2016_se <- sqrt(diag(law_obeyance_sport_active_all_2016_se))

#GR5
law_obeyance_sport_passive_active_2016<-lm(law_obeyance~member_sport_passive_active+
                                      age_group +income +edu+
                                      factor(region_name),
                                    data=data_2016, na.action=na.exclude)
summary(law_obeyance_sport_passive_active_2016)
law_obeyance_sport_passive_active_2016_se <- vcovHC(law_obeyance_sport_passive_active_2016, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_passive_active_2016_se <- sqrt(diag(law_obeyance_sport_passive_active_2016_se))

#GR6
law_obeyance_sport_passive_nonmb_2016<-lm(law_obeyance~member_sport_passive_nonmb+
                                            age_group +income +edu+
                                            factor(region_name),
                                          data=data_2016, na.action=na.exclude)
summary(law_obeyance_sport_passive_nonmb_2016)
law_obeyance_sport_passive_nonmb_2016_se <- vcovHC(law_obeyance_sport_passive_nonmb_2016, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_passive_nonmb_2016_se <- sqrt(diag(law_obeyance_sport_passive_nonmb_2016_se))

#GR7
law_obeyance_sport_passive_all_2016<-lm(law_obeyance~member_sport_passive_all+
                                          age_group +income +edu+
                                          factor(region_name),
                                        data=data_2016, na.action=na.exclude)
summary(law_obeyance_sport_passive_all_2016)
law_obeyance_sport_passive_all_2016_se <- vcovHC(law_obeyance_sport_passive_all_2016, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_passive_all_2016_se <- sqrt(diag(law_obeyance_sport_passive_all_2016_se))

#GR1
law_obeyance_sport_2010<-lm(law_obeyance~member_sport+
                              age_group +income +edu+
                              factor(region_name),
                            data=data_2010, na.action=na.exclude)
summary(law_obeyance_sport_2010)
law_obeyance_sport_2010_se <- vcovHC(law_obeyance_sport_2010, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_2010_se <- sqrt(diag(law_obeyance_sport_2010_se))

#GR2
law_obeyance_sport_active_passive_2010<-lm(law_obeyance~member_sport_active_passive+
                                             age_group +income +edu+
                                             factor(region_name),
                                           data=data_2010, na.action=na.exclude)
summary(law_obeyance_sport_active_passive_2010)
law_obeyance_sport_active_passive_2010_se <- vcovHC(law_obeyance_sport_active_passive_2010, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_active_passive_2010_se <- sqrt(diag(law_obeyance_sport_active_passive_2010_se))

#GR3
law_obeyance_sport_active_nomb_2010<-lm(law_obeyance~member_sport_active_nomb+
                                          age_group +income +edu+
                                          factor(region_name),
                                        data=data_2010, na.action=na.exclude)
summary(law_obeyance_sport_active_nomb_2010)
law_obeyance_sport_active_nomb_2010_se <- vcovHC(law_obeyance_sport_active_nomb_2010, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_active_nomb_2010_se <- sqrt(diag(law_obeyance_sport_active_nomb_2010_se))

#GR4
law_obeyance_sport_active_all_2010<-lm(law_obeyance~member_sport_active_all+
                                         age_group +income +edu+
                                         factor(region_name),
                                       data=data_2010, na.action=na.exclude)
summary(law_obeyance_sport_active_all_2010)
law_obeyance_sport_active_all_2010_se <- vcovHC(law_obeyance_sport_active_all_2010, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_active_all_2010_se <- sqrt(diag(law_obeyance_sport_active_all_2010_se))

#GR5
law_obeyance_sport_passive_active_2010<-lm(law_obeyance~member_sport_passive_active+
                                      age_group +income +edu+
                                      factor(region_name),
                                    data=data_2010, na.action=na.exclude)
summary(law_obeyance_sport_passive_active_2010)
law_obeyance_sport_passive_active_2010_se <- vcovHC(law_obeyance_sport_passive_active_2010, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_passive_active_2010_se <- sqrt(diag(law_obeyance_sport_passive_active_2010_se))

#GR6
law_obeyance_sport_passive_nonmb_2010<-lm(law_obeyance~member_sport_passive_nonmb+
                                            age_group +income +edu+
                                            factor(region_name),
                                          data=data_2010, na.action=na.exclude)
summary(law_obeyance_sport_passive_nonmb_2010)
law_obeyance_sport_passive_nonmb_2010_se <- vcovHC(law_obeyance_sport_passive_nonmb_2010, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_passive_nonmb_2010_se <- sqrt(diag(law_obeyance_sport_passive_nonmb_2010_se))

#GR7
law_obeyance_sport_passive_all_2010<-lm(law_obeyance~member_sport_passive_all+
                                          age_group +income +edu+
                                          factor(region_name),
                                        data=data_2010, na.action=na.exclude)
summary(law_obeyance_sport_passive_all_2010)
law_obeyance_sport_passive_all_2010_se <- vcovHC(law_obeyance_sport_passive_all_2010, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_passive_all_2010_se <- sqrt(diag(law_obeyance_sport_passive_all_2010_se))

#GR1
law_obeyance_sport<-lm(law_obeyance~member_sport+
                         age_group +income +edu+
                         factor(region_name),
                       data=combined, na.action=na.exclude)
summary(law_obeyance_sport)
law_obeyance_sport_se <- vcovHC(law_obeyance_sport, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_se <- sqrt(diag(law_obeyance_sport_se))

#GR2
law_obeyance_sport_active_passive<-lm(law_obeyance~member_sport_active_passive+
                                        age_group +income +edu+
                                        factor(region_name),
                                      data=combined, na.action=na.exclude)
summary(law_obeyance_sport_active_passive)
law_obeyance_sport_active_passive_se <- vcovHC(law_obeyance_sport_active_passive, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_active_passive_se <- sqrt(diag(law_obeyance_sport_active_passive_se))

#GR3
law_obeyance_sport_active_nomb<-lm(law_obeyance~member_sport_active_nomb+
                                     age_group +income +edu+
                                     factor(region_name),
                                   data=combined, na.action=na.exclude)
summary(law_obeyance_sport_active_nomb)
law_obeyance_sport_active_nomb_se <- vcovHC(law_obeyance_sport_active_nomb, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_active_nomb_se <- sqrt(diag(law_obeyance_sport_active_nomb_se))

#GR4
law_obeyance_sport_active_all<-lm(law_obeyance~member_sport_active_all+
                                    age_group +income +edu+
                                    factor(region_name),
                                  data=combined, na.action=na.exclude)
summary(law_obeyance_sport_active_all)
law_obeyance_sport_active_all_se <- vcovHC(law_obeyance_sport_active_all, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_active_all_se <- sqrt(diag(law_obeyance_sport_active_all_se))

#GR5
law_obeyance_sport_passive_active<-lm(law_obeyance~member_sport_passive_active+
                                 age_group +income +edu+
                                 factor(region_name),
                               data=combined, na.action=na.exclude)
summary(law_obeyance_sport_passive_active)
law_obeyance_sport_passive_active_se <- vcovHC(law_obeyance_sport_passive_active, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_passive_active_se <- sqrt(diag(law_obeyance_sport_passive_active_se))

#GR6
law_obeyance_sport_passive_nonmb<-lm(law_obeyance~member_sport_passive_nonmb+
                                       age_group +income +edu+
                                       factor(region_name),
                                     data=combined, na.action=na.exclude)
summary(law_obeyance_sport_passive_nonmb)
law_obeyance_sport_passive_nonmb_se <- vcovHC(law_obeyance_sport_passive_nonmb, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_passive_nonmb_se <- sqrt(diag(law_obeyance_sport_passive_nonmb_se))

#GR7
law_obeyance_sport_passive_all<-lm(law_obeyance~member_sport_passive_all+
                                     age_group +income +edu+
                                     factor(region_name),
                                   data=combined, na.action=na.exclude)
summary(law_obeyance_sport_passive_all)
law_obeyance_sport_passive_all_se <- vcovHC(law_obeyance_sport_passive_all, type = "HC0", cluster = "PSU_name")
law_obeyance_sport_passive_all_se <- sqrt(diag(law_obeyance_sport_passive_all_se))


#######
#FINAL#
#######


return_graphs_coef<-function(models_run, models_run_alias){
  #Process:
  #-takes list of regression objects, associates them with a 
  #list of "official" label, keep the "treat" variable and
  #produce ggplots with those
  
  #Input:
  #-List of Regression model object after lm
  #-List of label for the DV that you want displayed in ggplot
  
  #Output:
  #-A ggplot with regression coefficients and robust CI
  
  names(models_run)<-models_run_alias
  df<-list()
  length(df)<-length(models_run)
  counter = 1
  for (i in 1:length(models_run)){
    #print(i)
    x<-models_run[[i]]
    name_x<-names(models_run)[i]
    #print(name_x)
    model_tidy<-rerun_regression(x, name_x)
    model_tidy$counter<-counter
    #print(model_tidy)
    df[[i]] = model_tidy
    counter = counter +1}
  print(df)
  big_data = do.call(rbind, df)
  irislabs1 <- rev(big_data$dv_official)
  irislabs2 <- rev(big_data$no_obs)
  #return(big_data)
  mean_feature <- ggplot(data = big_data, aes(x=estimate, y=as.numeric(reorder(dv_official, desc(counter))))) +
    geom_point(size = 2) +
    geom_errorbar(aes(xmin = lb_95, xmax = ub_95), size = 0.9, width = 0)+
    scale_y_continuous(name = "",
                       breaks = 1:length(irislabs1),
                       labels = irislabs1,
                       sec.axis = sec_axis(~.,
                                           breaks = 1:length(irislabs2),
                                           labels = irislabs2))+
    geom_vline(xintercept = 0) +
    theme_bw() +
    #scale_y_discrete(name = "") +
    scale_x_continuous(name = "Estimate") +
    #ggtitle(glue("{title_spec}")) +
    theme(plot.title = element_text(hjust = 0.5),
          axis.text.x = element_text(size=14),
          axis.text.y = element_text(size=14),
          axis.title=element_text(size=16))  
  #return(big_data)
  return(mean_feature)
}


rerun_regression<-function(regression_name, string_dv){
  #Process:
  #-takes already existing regression objects, associates them with an 
  #"official" label, calculates robust SE and CI and puts them in a tibble
  
  #Input:
  #-Regression model object after lm
  #-The label for the DV that you want displayed in ggplot
  
  #Output:
  #-A tibble with regression coefficients from the model that you ran with:
  #--robust SE
  #--robust CI
  
  #Recreate dataframe from the regression object
  remake_data<-regression_name$model
  #Scale the DV
  remake_data$dv_scaled<-scale(remake_data[1])
  #Take the string of the DV and IV
  dv<-regression_name$terms[[2]]
  ivs<-regression_name$terms[[3]]
  ivs_str<-paste(ivs)[2]
  dv_str<-paste(dv, "~", collapse = "")
  dv_ivs_str<-paste(dv_str, ivs_str, sep=" ")
  dv_str_scaled<-paste("dv_scaled", "~", collapse = "")
  dvscales_ivs_str<-paste(dv_str_scaled, ivs_str, sep=" ")
  #Create specification and run model
  specification = as.formula(dvscales_ivs_str)
  model <- lm(specification, data = remake_data)
  no_obs<-nobs(model)
  no_obs<-paste("Obs: ", as.character(no_obs), sep = "")
  #Calculate robust SE
  model_cov <- vcovHC(model, type = "HC0", cluster = "PSU_name")
  model_se <- sqrt(diag(model_cov))
  #Create tibble
  model_tidy<-tidy(model)
  model_tidy$no_obs<-no_obs
  
  model_se_tidy<-tidy(model_se)
  names(model_se_tidy)[1] <- "term"
  names(model_se_tidy)[2] <- "robust.se"
  model_tidy2<-left_join(model_tidy, model_se_tidy, by = c("term"="term"))
  
  #Calculate CI
  model_tidy3<-model_tidy2
  model_tidy3$lb_95<-model_tidy3$estimate-1.96*model_tidy3$robust.se
  model_tidy3$ub_95<-model_tidy3$estimate+1.96*model_tidy3$robust.se
  
  #Only keep the treat variable
  #model_tidy4<-subset(model_tidy3, term=="treat")
  model_tidy4<-model_tidy3[model_tidy3$term %like% "member", ]
  model_tidy4$dv_official<-string_dv
  return(model_tidy4)
}


clean_col_labels<-function(column_labels){
  #Process:
  #-takes list of column labels for stargazer
  #and removes the unnecessary tags
  
  #Input:
  #-list of column labels for stargazer
  
  #Output:
  #-retruns list of clean labels
  df<-list()
  length(df)<-length(column_labels)
  counter = 1
  for (i in column_labels){
    #print(i)
    j<-str_remove(i, '\\{')
    k<-str_remove(j, '\\}')
    l<-str_remove(k, '\\\\')
    m<-str_remove(l, 'shortstack')
    n<-gsub("\\\\", "\n", m)
    o<-gsub("\n\n", "\n", n)
    p<-gsub(",", "", o)
    df[counter] = p
    counter = counter +1
    df<-unlist(df)
  }
  return(df)
}



#############
#Final Table#
#############

#Active and Passive vs. Non-Members
#Active vs. Passive
#Passive vs. Non-members

#law_obeyance
models_run<-list(law_obeyance_sport_2016,
                 law_obeyance_sport_2010,
                 law_obeyance_sport_active_passive_2016,
                 law_obeyance_sport_active_passive_2010,
                 law_obeyance_sport_passive_nonmb_2016,
                 law_obeyance_sport_passive_nonmb_2010)


column_labels<-c("\\shortstack{Active and Passive\\\\ vs. Non-Members}",
                 "",
                 "\\shortstack{Active\\\\ vs. Passive}",
                 "",
                 "\\shortstack{Passive\\\\ vs. Non-Members}",
                 "")

models_run_alias_clean<-clean_col_labels(column_labels)



names(models_run)<-models_run_alias_clean
df<-list()

length(df)<-length(models_run)
counter = 1
for (i in 1:length(models_run)){
  print(i)
  x<-models_run[[i]]
  name_x<-names(models_run)[i]
  #print(name_x)
  model_tidy<-rerun_regression(x, name_x)
  model_tidy$counter<-counter
  model_tidy$year<-ifelse(stri_detect_fixed(as.character(models_run[[i]]$call$data), "2010"), "2010", "2016")
  print(model_tidy)
  df[[i]] = model_tidy
  counter = counter +1}

big_data = do.call(rbind, df)
irislabs1 <- rev(big_data$dv_official)
irislabs2 <- rev(big_data$no_obs)
big_data$year<-factor(big_data$year, levels=c("2010", "2016"))



mean_feature <- ggplot(data = big_data, aes(x=estimate, 
                                            y=rev(counter),
                                            colour= year))+
  geom_point(aes(shape=year), position = position_dodge(0.5), size = 4) +
  geom_errorbar(aes(xmin = lb_95, xmax = ub_95), size = 0.9, width = 0)+
  scale_y_continuous(name = "",
                     breaks = 1:length(irislabs1),
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks = 1:length(irislabs2),
                                         labels = irislabs2))+
  scale_colour_manual(name = "Sample", 
                      label = c("2010", "2016"),
                      values=c("grey50", "grey5"))+
  scale_shape_manual(name = "Sample", 
                     label = c("2010", "2016"),
                     values=c(15, 16))+
  geom_vline(xintercept = 0) +
  theme_bw() +
  #scale_y_discrete(name = "") +
  scale_x_continuous(name = "Estimate", breaks=seq(-0.2,0.6,0.2), limits = c(-0.2,0.6))+
  #ggtitle(glue("{title_spec}")) +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title=element_text(size=16),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))

mean_feature2<-reposition_legend(mean_feature, 'top right')



ggsave(mean_feature2, file = './paper/graphs/figure_2.jpg', 
       height=16.42, width=18.2, units = "cm", dpi=300)

